Numerical study of the phase transitions in the 
two-dimensional Z(5) vector model 



0. Borisenko* 

Bogolyubov Institute for Theoretical Physics, 
National Academy of Sciences of Ukraine, 
03680 Kiev, Ukraine 

G. Cortese**, R. Fiore** 

Dipartimento di Fisica, Universita della Calabria, 
and Istituto Nazionale di Fisica Nucleare, Gruppo collegato di Cosenza 
1-87036 Arcavacata di Rende, Cosenza, Italy 

M. Gravina*** 

Laboratoire de Physique Theorique, Universite de Paris-Sud 11, Batiment 210 

91405 Orsay Cedex France 
and Department of Physics, University of Cyprus, P.O. Box 20357, Nicosia, Cyprus 

A. Papa** 

Dipartimento di Fisica, Universita della Calabria, 
and Istituto Nazionale di Fisica Nucleare, Gruppo collegato di Cosenza 
1-87036 Arcavacata di Rende, Cosenza, Italy 

Abstract 

We investigate the critical properties of the two-dimensional Z(5) vector model. 
For this purpose, we propose a new cluster algorithm, valid for Z(N) models with 
odd values of N. The two-dimensional Z(5) vector model is conjectured to exhibit 
two phase transitions with a massless intermediate phase. We locate the position of 
the critical points and study the critical behavior across both phase transitions in 
details. In particular, we determine various critical indices and compare the results 
with analytical predictions. 
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1 Introduction 



The Berezinskii-Kosterlitz-Thouless (BKT) phase transition is known to take place in a 
variety of two-dimensional (2D) systems: certain spin models, two-dimensional Coulomb 
gas, sine-Gordon model, Solid-on- Solid model, etc., the most popular and elaborated case 
being the two-dimensional XY model [H El E] . There are several indications that this type 
of phase transition is not a rare phenomenon in gauge models at finite temperature: one 
can argue that in some three-dimensional lattice gauge models the deconfinement phase 
transition is of BKT type as well. Here we are going to study an example of lattice spin 
model where this type of transition exhibits itself, namely the 2D Z(N) spin model, also 
known as vector Potts model. 

Consider a 2D lattice A = L 2 with linear extension L and impose periodic boundary 
conditions on spin fields in both directions. The partition function of the model can be 
written as 



N-1 



zgA s(x)=0 



n n q - s ( x + e ™)) 



.xeA 71=1,2 



(1) 



In the standard formulation the most general Z(AQ-invariant Boltzmann weight with 
N — 1 different couplings is 



Q(s) 



exp 



N-1 , 

> t Pk cos — —s 

.k=l 



N 



(2) 



In the Villain formulation the Boltzmann weight reads instead 



QM 



E 



exp 
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Some details of the critical behavior of 2D Z(N) spin models are well known - see 
the review in Ref. jl]. The Z(N) spin model in the Villain formulation has been 
studied analytically in Refs. EJ [71 El [9]. It was shown that the model has at least 
two phase transitions when N > 5. The intermediate phase is a massless phase with 
power-like decay of the correlation function. The critical index rj has been estimated 
both from the renormalization group (RG) approach of the Kosterlitz-Thouless type and 
from the weak-coupling series for the susceptibility. It turns out that r]^^) = 1/4 at 
the transition point from the strong coupling (high-temperature) phase to the massless 

(2) 

phase, i.e. the behavior is similar to that of the XY model. At the transition point /3c 
from the massless phase to the ordered low-temperature phase one has r]{(5 { c 2) ) = A/N 2 . 
A rigorous proof that the BKT phase transition does take place, and so that the massless 
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phase exists, has been constructed in Ref. [TD] for both Villain and standard formulations 
(with one non- vanishing coupling 0x). Monte Carlo simulations of the standard version 
with N = 6,8, 12 were performed in Ref. [11] . Results for the critical index rj agree well 
with the analytical predictions obtained from the Villain formulation of the model. 

In this paper we thoroughly investigate the case N = 5, the lowest number where the 
BKT transition is expected. Precisely, we concentrate on the standard formulation (J2J) 
with one non-zero coupling The motivation of our study is three- fold: 

1. to compute critical indices at the transition points, which could serve as checking 
point of universality; 

2. to shed light on the discrepancy in the literature concerning the Z(5) model; 

3. to develop and test a new version of Monte Carlo cluster algorithm valid for odd 
values of N. 

The first motivation is related to the study of the finite-temperature transitions in 
3D Z(N) and SU(N) lattice gauge theory (LGT). It is expected that in 3D Z(N) LGT 
a deconfinement phase transition takes place at finite temperature. There is no precise 
statement about the order of the phase transition, but presumably it is of the BKT type 
if N > 4. If it is the case, the Svetitsky-Yaffe conjecture [12] implies that the 3D Z(N) 
LGT is in the universality class of 2D vector Potts model. Moreover, it can be proven 
that, in the strong coupling region with respect to the spatial coupling, the 3D Z(N) 
LGT reduces to a 2D Z(N) model with the general Boltzmann weight (jSJ), and such that 
/?i = /?4 /?2 = /^3, for N = 5. Here, (3k are effective couplings which depend on the 
gauge coupling and the temporal extension N t . Thus, our Z(5) model represents a good 
approximation to 3D Z(5) LGT in this region. 

Next, let W(x) e SU(N) and consider the following effective action in 2D 



The effective action (jl]) can be regarded as the simplest effective model for the Polyakov 
loop which can be derived in the strong coupling region of 3D SU(N) LGT at finite 
temperature. It possesses Z(N) global symmetry and thus may well exhibit the BKT 
transitions which belong to the universality class of the corresponding vector Potts model. 
Therefore, our investigation here can be viewed as a preliminary step in studying decon- 
finement phase transition in 3D Z(N) and SU(N) LGTs. 

The second motivation reflects the fact that many features of the critical behavior of 
the Z(5) model are not reliably established. Moreover, there are certain discrepancies 
even in determining the nature of the phase transition, i. e. whether the phase transition 
is of BKT type or not. 




(4) 
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Let us briefly summarize the present state of affairs. 

• The rigorous proof of the massless phase existence in 2D Z(N) models utilizes 
methods which do not allow to establish the exact value of N above which the BKT 
phase transition exists [TO"] . 

• Reliable analytical calculations can only be performed with the Villain formulation 
([2D; the RG study of Ref. [S] predicts that the massless phase exists for all iV > 4. 

• Some information on the phase structure of the general Z(N) spin models can be 
obtained through the duality transformations (see, e.g. p3]). These transformations 
cannot be used to establish the position of the critical points in the Z(b) model [2]. 
However, duality transformations relate the two critical points and thus can be used 
to verify the accuracy of numerical data. Moreover, one can predict an approximate 
phase diagram and argue that the massless phase and the BKT transition exist for 
N = 5 in a certain region of the parameter space [XSJ [16] . Important in this context 
is the rigorous proof of the existence of the massless phase, constructed in Ref. [15J 
for the generalized Villain formulation. This generalized formulation contains the 
vector Potts model defined in (JSJ) with one non-zero coupling j3\ as a particular case. 

For completeness we mention that in Refs. [I7J it was suggested that there is only one 
first order phase transition in Z(5) model. This, however, contradicts the rigorous 
results of [T5] . 

• An analytical prediction for the critical index r] has been obtained for the Villain 
formulation in Ref. [5]: ^(A^) = 1/4 and rj{fi^) = 0.16 (for N = 5). The situation 
remains unclear for the index v which governs the behavior of the correlation length. 
Normally, the value of v can be estimated from the solution of the RG equations 
like in the XY model [3J. The analytical solution of the system of RG equations for 
Z(N) vector models is unknown (see Ref. [5]). Therefore, strictly speaking, there 
are no strong theoretical arguments indicating that v = 1/2, similarly to the XY 
model. Moreover, the strong coupling expansion of the Z(5) model combined with 
Pade approximants predicts that v fa 0.22 [5]. 

• Monte Carlo simulations of the Z(N) model have been performed in Refs. [TTJ [18j 
\19\ [20] . Results of Ref. [18] , though obtained on rather small lattices, indicate that 
the BKT transition takes place in models with N > 8. This contradicts the results 
of [HI [20] which well agree with the BKT behavior for N = 6. 

However, most recent simulations of the helicity modulus in the Z(h) model at 0c 
do not agree with what is expected at the BKT transition [TjJ]. Namely, at the 
critical point the helicity modulus is expected to jump discontinuously to zero, and 
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the jump is observed for Z(6) model, while in Z(5) the helicity modulus stays small 
but non- vanishing in the high-temperature region f3 < /3c . 

Still, it remains unclear how this behavior of the helicity modulus influences other 
features of the BKT transition. The key feature of the massless BKT phase in Z(N) 
models is the enhancement of the discrete symmetry of the Hamiltonian: the symmetry 
of the ground state in the intermediate phase is rather U(l) than Z(N) [TO]. This can be 
seen in a characteristic distribution of the complex magnetization, in the power-like decay 
of the correlation functions in the massless phase, in the vanishing of the beta-function, 
etc.. Also, on the basis of the universality one could conjecture that the critical indices in 
the Z{5) model are the same both in the standard and Villain formulations. We are not 
aware of any numerical calculations of these quantities in Z(5). Here we would like to fill 
this gap by computing various quantities and extracting critical indices at both transitions. 
Among the main quantities calculated in this paper there are Binder cumulants. As RG 
invariant quantities, Binder cumulants are very useful in locating the critical couplings 
and determining the nature of the phase transition. Indeed, the computation of Binder 
cumulants proved to be very efficient in studying BKT transitions in a variety of models, 
like the XY model, the discrete Gaussian model and the SOS model [2T] . We therefore 
believe that these cumulants are of great value also in the investigation of phase transitions 
in 2D Z(N) models. Preliminary results of our study have been presented in Ref. [22J. 

The paper is organized as follows: in Section [5] we describe the set-up of the Monte 
Carlo simulation and the newly developed cluster algorithm; in Section 3 we introduce the 
observables adopted in this work and study the transition from the high-temperature to 
the massless phase; in Section 4 we move on to consider the transition from the massless 
to the low-temperature ordered phase; finally in Section 5 we draw our conclusions. In 
the Appendix we check the consistency of our determination of the critical couplings with 
the duality transformations. 

2 Algorithm and numerical set-up 

In this work we concentrate our attention to the model defined by Eqs. ([T]) and 02]), with 
only one non-zero coupling, 0i = /3 Q. This model is known in the literature also as 
./V-state ferromagnetic clock model and is a discrete version of the continuous XY (plane 
rotator) model. It consists of 2D planar spins restricted to N evenly spaced directions, 
with spin interaction energy proportional to their scalar product. 

1 All forthcoming tables and plots refer to the case N = 5, but we nevertheless present all definitions 
and formulae for a generic N. 
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The Hamiltonian of the model is 

# = cos , Si = 0,l,...,N-l, (5) 

with summation taken over nearest-neighbor sites. For N = 2 this is the Ising model, 
whereas in the N — > oo limit we get the XY model. 

Here we develop a new algorithm, valid for odd N, by which an accurate numerical 
study of the model can be performed for N = 5, i.e. the smallest N value for which the 
phase structure described in the Introduction holds. 

Here are the steps of our cluster algorithm for the update of a spin configuration {si}: 

• choose randomly n in the set {0, 1, 2, . . . , N — 1} 

• build a cluster configuration according to the following probability of bond activation 
between neighboring sites ij 

f 1 — exp(— 28 aiCtj) if a^aj > ,,. . f 2n \ 

& ={ otherwise ' Wlth " fc = Sm { N {Sk ~ n) ) 

• "flip" each cluster, with probability 1/2, by replacing all its spins according to the 
transformation 

Si ->■ mod(-Sj + 2n + N, N) , 

which amounts to replacing each spin in a cluster by the spin Sj for which otj = 
—oii\ equivalently, if the spins are mapped into the N roots of unity in the complex 
plane, the above replacement means flipping the component of each spin transverse 
to the direction identified by n. 

It is easy to prove that this cluster algorithm fulfills the detailed balance. 

We have tested the efficiency of the cluster algorithm against the standard heat-bath 
algorithm. On a lattice with L = 64 we simulated the model with N = 5 and deter- 
mined the autocorrelation time r of three observables: the energy, defined in Eq. ((Sj), 
the magnetization M L and the population Sl, to be defined below. We considered three 
B values (0.80, 1.10 and 1.50) lying, respectively, in the high-temperature, massless and 
low-temperature phase of the model. Results are summarized in Table [TJ whose last two 
columns give also the number of sweeps needed to reach thermal equilibrium and the 
computer time to collect a 50k statistics. 

At /3=0.80 and 1.10 the autocorrelation time in the cluster algorithm is lower than 
in the heat-bath for the energy and much lower for magnetization and population. At 
8 = 1.50, deep in the low-temperature ordered phase, r is systematically higher in the 
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Table 1: Cluster versus heat-bath in Z(5) on a 64 2 lattice at three values of /3: auto- 
correlation time r for three observables (energy, magnetization Ml and population Sl), 
number of thermalization sweeps and computer time for 50k updates. 





Energy Ml Sl 


Thermalization Time 


(3 = 0.80 


cluster 5.135(82) 1.528(38) 1.494(39) 
heat-bath 5.43(33) 12.83(27) 12.65(28) 


~ 10 46.69 s 
~ 100 1290 s 


p = 1.10 


cluster 7.36(11) 5.56(10) 7.18(11) 
heat-bath 10.11(48) 48.3(4.8) 60.6(6.1) 


~ 10 45.14 s 
~ 1000 1194 s 


(3 = 1.50 


cluster 8.97(17) 8.71(17) 8.84(17) 
heat-bath 2.38(12) 3.73(15) 3.78(13) 


~ 100 42.60 s 
~ 6500 1064 s 



P = 0.80 
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Re(M ) 
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-0.6 
-0.8 
-1 



P = 1.10 



1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 
Re(M ) 




-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 
Re(M ) 



Figure 1: (Color online) Scatter plot of the complex magnetization Ml at (3 = 
0.80, 1.10, 1.50 in Z{h) on a 64 2 lattice. 



cluster than in the heat-bath. This is a consequence of the lowering of the bond activation 
probability for increasing f3. This drawback, however, is compensated by the higher simu- 
lation speed, with respect to the heat-bath algorithm. Moreover since the two transitions 
in the 2D Z(5) model are rather close (see below), there is no doubt that the cluster 
algorithm is strongly preferable. 

The improvement brought along by the cluster algorithm becomes more visible when 
the dynamical critical exponent z is considered, defined as r ~ £ z , where £ is the 
correlation length. We have evaluated z in the 2D Z(5) model on lattices with L = 
16, 32, 64, 128, 256, 384, 512 at both transition points, using the autocorrelation time of 
the magnetization Ml- Since at both points the correlation length diverges, the expected 
scaling law becomes r ~ L z . We got in all cases that r keeps almost constant at w 7, 
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thus implying z ~ 0, i.e. no critical slowing down. 

The three phases exhibited by the 2D Z(5) spin model can be characterized by means 
of two observables: the complex magnetization Ml and the population Sl- 

The complex magnetization is given by 

Ml = 7^E ex P Oji 3 ) = W e# ■ ( 6 ) 

In Fig. [T]we show the scatter plot of Ml on a lattice with L = 64 in Z(5) at three values of 
0, each representative of a different phase: = 0.80 (high-temperature, disordered phase), 
= 1.10 (BKT massless phase) and = 1.50 (low-temperature, ordered phase). As we 
can see we pass from a uniform distribution (low 0) to a ring distribution (intermediate 
0) and finally to five isolated spots (high 0). 

The naive average of the complex magnetization gives constantly zero, therefore Ml 
is not an order parameter. An observable to detect the transition from one phase to the 
other is instead the absolute value \Ml\ of the complex magnetization. In Fig. [2] we show 
the behavior of \M L \ and of its susceptibility, 

X [ M) = L 2 ((|M L | 2 ) - (|M L |) 2 ) , (7) 

in Z(5) on lattices with L ranging from 16 to 1024 over a wide interval of values. On 
each lattice the susceptibility Xl clearly exhibits two peaks, the first of them, more 
pronounced than the other, identifies the pseudocritical coupling 0pc(L) at which the 
transition from the disordered to the massless phase occurs, whereas the second corre- 
sponds to the pseudocritical coupling 0pJ(L) of the transition from the massless to the 
ordered phase. It is evident from Fig. [2] that \Ml\ is particularly sensitive to the first 
transition, thus making this observable the best candidate for studying its properties. 

As order parameter to better detect the second transition, i. e. that from the massless 
to the ordered phase, we chose instead the population Sl, defined as 



a v 



N 



max i= o,jv_i(7Xi) 



L 2 N 



(8) 



where represents the number of spins of a given configuration which are in the state 
Si. In a phase in which there is not a preferred spin direction in the system (disorder), 
we have n« ~ L 2 /N for each index i, therefore Sl ~ 0. Otherwise, in a phase in which 
there is a preferred spin direction (order), we have Ui ~ L 2 for a given index i, therefore 
S l ~ 1- In Fig. [3] we show the behavior of Sl and of its susceptibility 

= L 2 ((Sl) - (S L ) 2 ) , (9) 



S 




p p 



Figure 2: (Color online) Behavior of \Ml\ (left) and of its susceptibility Xl (right) versus 
j3 in Z(5) on lattices with L ranging from 16 to 1024. 



in Z(5) on lattices with L ranging from 16 to 1024 over a wide interval of (5 values. Again 
the peaks signalling the two transitions are clearly visible and their positions agree with 
Fig. El but now the second one is more pronounced. 

Other observables which have been used in this work are the following: 

• the real part of the "rotated" magnetization, Mr = \Ml\ cos(Nip), 

• the order parameter introduced in Ref. [20], = cos(Nijj), 

where ip is the phase of the complex magnetization defined in Eq. ([6]). 

In the next two Sections we will study separately the two transitions of 2D Z(5) and 
determine some of the related critical indices. For all observables considered in this work 
we collected typically 100k measurements, on configurations separated by 10 updating 
sweeps. For each new run the first 10k configurations were discarded to ensure ther- 
malization. Data analysis was performed by the jackknife method over bins at different 
blocking levels. 

3 The transition from the high-temperature to the 
massless phase 

The first inflection point in the plot of the magnetization |Ml| and the first peak in the 
plot of the susceptibility Xl ( see Fig- EJ) indicate the transition from the disordered to the 
massless phase. The couplings where this transition occurs (denoted as the pseudocritical 
couplings j3pc(L)) have been determined by a Lorentzian interpolation around the peak 
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Figure 3: (Color online) Behavior of Sl (left) and of its susceptibility xf (right) versus 
/3 in Z(5) on lattices with L ranging from 16 to 1024. 



of the susceptibility Xl^ • Their values are summarized in the second column of Table [2j 
We observe that, when the lattice size L grows, (3pc(L) increases towards the infinite 
volume critical coupling /3c and that the susceptibility Xl S oes ^° zero ^ ess ra pidly for 
/3 > (3pc , as expected in the BKT scenario. 

In order to apply the finite size scaling (FSS) program, the location of the infinite 
volume critical coupling (3^ is needed. In Refs. [231 121] this was done by extrapolating 
the pseudocritical couplings to the infinite volume limit, according to a suitable scaling 
law. First order transition is ruled out by data in Table [2] Second order transition, 
though not incompatible with data in Table [2J is to be excluded, due to the vanishing 
of the long distance correlations combined with the clusterization property (we will come 
back to this point in the last Section). Therefore, we assume that the transition is of BKT 
type and adopt the scaling law dictated by the essential scaling of the BKT transition, 
i.e. £ ~ e bt , which reads 



PS? = 0c W + r 

P (lnL + 5)^ 



v 



(10) 



1/2 



The index v characterizes the universality class of the system. For example, 
holds for the 2D XY universality class. 

Unfortunately, 4-parameter fits of the data for /3p C \L) give very unstable results for 
the parameters. This led us to move to 3-parameter fits of the data, with v fixed at 1/2. 
We found, as best fit with the MINUIT optimization code, 



(3^ = 1.0602(20) , A x = -2.09(20) , B x = 0.27(18) , x7 d -°- f - = °- 48 > L " 



64 
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Table 2: Values of (3$ in Z(5) on L 2 lattices. The last two columns give the susceptibility 
Xl and the magnetization \Ml\ at the infinite volume coupling constant /3<T =1.0510. 



L 






M L ($ i} ) 


16 


0.8523(20) 






32 


0.91429(90) 






64 


0.95373(40) 






128 


0.98054(30) 






256 


0.99838(20) 






384 


1.00621(10) 


187.9(1.2) 


0.48929(13) 


512 


1.01112(20) 


311.5(2.0) 


0.47181(13) 


640 




458.6(3.4) 


0.45918(11) 


768 




631.3(4.2) 


0.44863(11) 


896 




824.4(5.2) 


0.44004(11) 


1024 


1.01991(10) 


1040.0(6.9) 


0.43277(11) 



We observe that /3c is rather far from the value of (3$ on the largest available lattice, 
thus casting some doubts on the reliability of the extrapolation to the thermodynamic 
limit. For this reason, we turned to an independent method for the determination of /3c , 
based on the use of Binder cumulants. 

In particular, we considered the reduced 4-th order Binder cumulant defined as 

V {M) _ _ (\M L \ A ) 

and the cumulant B^ Mr ^ defined as 

p (M R) _ (\M R -(M R )\*) 

^ (\M R -{M R )\Y' 1 } 

Plots of the various Binder cumulants versus /3 show that data obtained on different lattice 
volumes align on curves that cross in two points, corresponding to the two transitions (see 
Figs. H]and[5]). We used also the reduced 4-th order Binder cumulant of the action which 
showed no crossing points nor volume-dependent dips, thus confirming the absence of first 
order phase transitions. 

We determined the crossing point by plotting the Binder cumulants versus (/3 — 
/3 c )(lnL) 1// ^, with v fixed at 1/2, and by looking for the optimal overlap of data from 
different lattices, by the \ 2 method (see Fig. [6] for an example of this kind of plots). As a 
result of this analysis we arrived at the following estimate: /3c = 1.0510(10). We observe 
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P 

Figure 4: (Color online) Reduced 4-th order Binder cumulant versus f3 on lattices 

with L ranging from 16 to 1024. 

that 0c is not compatible with the infinite volume extrapolation of the corresponding 
pseudocritical couplings, thus confirming our previous worries about the safety of the 
infinite volume extrapolation of (3$ . It should be noted, however, that a fit to /3^ c \l) 
with the law (fTUj) and with {3^ fixed at 1.0510 and v fixed at 1/2 gives a good x 2 /d.o.f., 
if only the three largest volumes are considered in the fit. 

We have also tested the strong coupling prediction of Ref. [5] suggesting v = 0.22. 
In particular, we have plotted the Binder cumulant versus (/3 — /3 c )(ln L) 1 /^, with v 
fixed now at the candidate value 0.22. We have seen that, varying (3 C on a wide interval, 
the overlap among curves from different lattices is always rather poor. 

Table 3: Results of the fit to the data of \Ml\(/3c) with the scaling law ([TBI on L 2 lattices 
with L > L m ; n . 





A 


P/v 


X 2 /d.o.f. 


384 


1.0299(21) 


0.12508(32) 


1.3 


512 


1.0294(32) 


0.12501(47) 


1.7 


640 


1.0371(49) 


0.12610(71) 


0.40 


768 


1.0305(89) 


0.1252(13) 


0.021 



We are now in the position to extract other critical indices and check therefore the 
hyperscaling relation. According to the standard FSS theory, in a LxL lattice at criticality 
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Figure 5: (Color online) Binder cumulant B\ versus (3 on lattices with L ranging from 
128 to 1024. 



Table 4: Results of the fit to the data of xi (Pc ) with the scaling law (fT4l) on L 2 lattices 
with L > L mm . 





A 


7 /z/ 


X 2 /d.o.f. 


384 


0.00586(30) 


1.7438(80) 


0.060 


512 


0.00602(48) 


1.740(12) 


0.018 


640 


0.00598(81) 


1.741(20) 


0.025 


768 


0.0062(14) 


1.735(34) 


0.0063 



the equilibrium magnetization \Mi\ should obey the relation \Ml\ ~ L @' u , for sufficiently 
large L 0. We performed a fit to the data of \M L \((3c^) (reported in the last column of 
Table |2]) on all lattices with size L not smaller than a given L m [ n according to the scaling 
law 

\M L \=AL-M" (13) 

and summarized our results in Table [3j 

The FSS behavior of the susceptibility defined in Eq. flTJ) is given by x^ ~ L 1 ^ , 
where 7/z/ = 2 — 77 and 77 is the magnetic critical index. We performed a fit to the data of 
Xl (/3c ) (reported in the third column of Table |2]) on all lattices with size L not smaller 

2 The symbol j3 here denotes a critical index and not, obviously, the coupling of the theory. In spite 
of this inconvenient notation, we are confident that no confusion will arise, since it will be always clear 
from the context which j3 is to be referred to. 
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Figure 6: (Color online) Reduced 4-th order Binder cumulant versus /3 c )(ln L) 1 ^ , 
for /3 C = 1.0510 and v = 1/2 on lattices with L ranging from 128 to 1024. 



than a given L min according to the scaling law 

xf ) = (14) 

and summarized our results in Table HI As we can see, for all values of L min considered, 
the value of the magnetic index El 77 W = 2 — j/u is compatible with 1/4. Note also that 
the hyperscaling relation 7/1/ + 2/3 /v = d, where d is the dimension of the system, is 
always satisfied within the statistical error. 

An independent determination of the magnetic index 77 can be achieved by the ap- 
proach developed in Ref. [23]: an effective 77 index is defined, through the spin-spin cor- 
relation function T(R), according to 

with Rq chosen equal to 10, as in Ref. [23]. This quantity is constructed in such a way 
that it exhibits a plateau in R if the correlator obeys the law 

T (R) * ^ , (16) 

valid in the BKT phase, (3 > (3^. In Fig. [7] we show the behavior of r]^(R) at the infinite 
volume critical coupling /3 C (1) = 1.0510 on lattices with L = 384,512,640,768,896,1024. 



3 The notation W in 7/ means "at the infinite volume critical coupling of the first transition". 
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It turns out that a plateau develops at small distances when L increases and that the 
extension of this plateau gets larger with L, consistently with the fact that finite volume 
effects are becoming less important. The plateau value of 77^ can be estimated at about 
0.24. We checked that this result is stable under variation of the parameter Rq. The 
discrepancy with the expected value of 1/4 can be explained by the imperfect localization 
of the critical point and/or by the effect of logarithmic corrections [251 126] that we were 
not able to include in our analysis. 



4 The transition from the massless to the low-tem- 
perature ordered phase 

The second inflection point in the plot of the population Sl and the second peak in the 
plot of the susceptibility \ L ( see Fig- El) indicate the transition from the massless to the 
ordered phase. The couplings where this transition occurs (denoted as the pseudocritical 

(2) 

couplings pp C (L)) have been determined by a Lorentzian interpolation around the peak 
of the susceptibility xf '■ Their values are summarized in the second column of Table 
Available results in the literature [51 [11] suggest that the correlation length diverges 
according to essential scaling scenario when the critical point is approached from above. 
Our aim is to check the validity of this statement and then extract relevant indices charac- 
terizing the system at this transition. Again, first order transition is ruled out by data in 
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Table 5: Values of f3 pc in Z(5) on L 2 lattices. The last two columns give the suscepti- 
bility Xi^ an d the rotated magnetization Mr at the infinite volume coupling constant 
=1.1048. 



L 


f& 






16 


1.1323(19) 






32 


1.1363(11) 






64 


1.13212(60) 






128 


1.12875(66) 






256 


1.12290(16) 






384 


1.12103(50) 


47116(77) 


0.1618(18) 


512 


1.11912(28) 


80057(139) 


0.1575(19) 


640 




120777(229) 


0.1557(20) 


768 




169358(298) 


0.1517(19) 


896 




224879(339) 


0.1502(16) 


1024 


1.11596(38) 


288151(532) 


0.1473(18) 



Tabled] (and by the aforementioned analysis of the Binder cumulant of the action), while 
second order is not. We assume that a BKT transition is at work here and, therefore, that 
pseudocritical couplings scale with L according to the law ( TTOl) . As before, 4-parameter 

(2) 

fits of the data for (3 pc (L) are unstable and we moved to 3-parameter fits of the data, 
with v fixed at 1/2, finding that the parameter B 2 turns out to be compatible with zero, 
so that, in fact, a 2-parameter fit works well: 

(3 (2) = 1.1042(12) , A 2 = 0.578(41) , B 2 = 0. , xVd-o.f. = 0.61 , L min = 128 . 

Now /3c is not far from the value of f3 pc on the largest available lattice, thus supporting 
the reliability of the extrapolation to the thermodynamic limit. 



Table 6: Results of the fit to the data of Xl r \Pc ) with the scaling law f[T41) on L 2 
lattices with L > L min . 



^min 


A 


l/v 


X 2 /d-o.f. 


384 


0.799(11) 


1. 


3459(21) 


0.23 


512 


0.791(17) 


1. 


3473(32) 


0.19 


640 


0.784(28) 


1. 


3487(53) 


0.22 


768 


0.793(50) 


1. 


3470(92) 


0.39 
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Table 7: Results of the fit to the data of Mr(/3c^) with the scaling law (I13p on L 2 lattices 
with L > L min . 



L<mm 


A 


0/v 


X 2 /d.o.f. 


384 


0.281(26) 


0.093(14) 


0.15 


512 


0.288(41) 


0.096(22) 


0.18 


640 


0.322(77) 


0.112(36) 


0.12 


768 


0.30(12) 


0.102(60) 


0.19 



In order to localize the critical coupling f3c , we looked for the crossing point at higher 
(3 of the Binder cumulant B^ Ir ^ defined in Eq. ffl2l and repeated the analysis based on 
the optimal overlap of data points when they are plotted against (ft — /3 c )(ln L) 1//y , with 
v fixed at 1/2. The same procedure was carried on using also the observable m^, which 
is itself an RG-invariant quantity and shares therefore the same properties of a Binder 
cumulant (see Fig. [8] for the behavior of versus on various lattices, which shows two 
crossing points, the one at higher (3 corresponding to the transition from the massless to 
the ordered phase). This analysis led to the result (3c = 1.1048(10), which agrees with 
the infinite volume extrapolation of the corresponding pseudocritical couplings. 

We can now determine the ratios of critical indices (3/v and j/v as we did in the pre- 
vious Section. It should be noted, however, that the population Sl and its susceptibility 
are not suitable observables for this purpose, being defined in a non-local manner and, 
thus, not directly related to the two-point correlator. We use, instead, the rotated magne- 
tization Mr and its susceptibility Xl an d compare their values at the infinite volume 

(2) ______ 

critical coupling (3c (see the last two columns of Table [5]) with the scaling laws (fT3"|) 
and (TT4"|) . respectively. Results for (3/v and 7/1/ are summarized in Tables |6]and[71 

All the values for rft*' = 2 — 7/1/ given in Table |6] are in agreement with the prediction 
4/iV 2 , which gives 0.16 for N = 5. The hyperscaling relation 7/1/ + 2/3/v = d is always 
satisfied, within the statistical error. 

The determination of the magnetic critical index based on the effective 77 index defined 
in ffl5|) is plagued, in this region of values of (3, by a sizeable dependence on the choice of 
the arbitrary parameter R . The shape of the curves for r] e s(R) and the way they depend 
on Rq suggest that here logarithmic corrections to the scaling could be at work. However, 
our data are not accurate enough to include them reliably in our fits. 

We conclude this Section by presenting several examples of a posteriori check of con- 
sistency of our determinations for r/ 1 ) and r]( 2 \ The basic idea is to build plots in which 
we correlate two RG-invariant quantities and to check that sequences of data points, cor- 
responding to different values of (3, fall on a universal curve, irrespective of the lattice size 
L [23 [28]. 
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Figure 8: (Color online) Behavior of m^p versus on lattices with L ranging from 128 to 
1024. 

The first example is the plot of the rescaled susceptibility ^ 1r ^ U 1-2 against the Binder 
cumulant b[ Mr \ One can see from Fig. 0(a) that for r] = 0.26 ~ rj^ data points from 
different lattices fall on the same curve in the lower branch, corresponding to values in 
the region of the first transition; for i] = 0.16 ~ r/ 2 \ on the contrary, data points from 
different lattices fall on the same curve in the upper branch, corresponding to values in 
the region of the second transition (see Fig. E](b)). 

Another example is provided by the plot of the rescaled magnetization M^L^ 2 against 
m^. For 7] = 0.16 ~ r]^ 2 \ again data points from different lattices fall on the same curve 
(see Fig. (TO]). 

5 Discussion and conclusions 

In this paper we have presented a wealth of numerical data aimed at shedding light on 
the phase structure of the 2D Z(5) vector model. By means of a Monte Carlo cluster 
updating algorithm, designed to work for Z(N) models with odd N, we have outlined a 
scenario compatible with the existence of three phases: disorder (small /3), massless or 
BKT (intermediate 0), order (large 0). We have determined 

• the critical points 0^ and 0^ in the infinite volume limit, by means of the FSS of 
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Figure 9: (Color online) Correlation between \l U 1 ^ 2 and the Binder cumulant B\ 
for (a) rj = 0.26 and (b) rj = 0.16 on lattices with L ranging from 128 to 1024. For 
r) = 0.26 (a) data from different lattices tend to fall on a universal curve in the lower 
branch, corresponding to (3 values in the region of the first transition. For 77 = 0.16 
(b) data from different lattices tend to fall on a universal curve in the upper branch, 
corresponding to (3 values in the region of the second transition. 



suitable Binder cumulant s; 

• the critical indices (3/v and 7/1/ at the two critical points, by means of the FSS of 
suitable definitions of the magnetization and of its susceptibility. 

The determination of (3c has been cross-checked with the infinite volume extrapo- 
lation of the (volume dependent) pseudocritical couplings of the second transition, as- 
suming essential scaling. We have found the following values of the critical couplings: 
fiP = 1.0510(10) and f3 { c 2) = 1.1048(10). As mentioned in the Introduction, values of 
the critical points are related by the duality transformations. In the Appendix we check, 
via the duality, the accuracy of our predictions and show that our determination of the 
critical couplings is in rather good agreement with it. 

The determination of the index r/ 1 ) = 2 — 7/1/ at the first transition has been cross- 
checked with the effective 77 index method. The values of the index 7/1/ at both critical 
points agree well with theoretical predictions obtained for the Villain formulation [5] 
thus supporting the conjecture that both standard and Villain formulations are in the 
same universality class. The behavior of the complex magnetization as well as the two- 
point correlation function strongly indicate that the intermediate phase is a massless 
phase whose symmetry is U(l). Moreover, using spin-wave-vortex approximation and 
conventional perturbation theory one can calculate the two-point correlation function 
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Figure 10: (Color online) Correlation between M^L^I 2 and for i] = 0.16 on lattices 
with L ranging from 128 to 1024. 



analytically and extract the perturbative beta-function in the intermediate phase. It turns 
out that the beta-function vanishes in this phase - a property which further supports the 
presence of the BKT transition and massless phase. 

For completeness, we have tested a scenario in which both transitions are second order. 
In this case, the infinite volume extrapolation of the pseudocritical couplings should obey 

tf c ' 2) = + 4- , 

whereas the overlap method of the Binder cumulants should work when data are plotted 
against [fi — /3 C )L». Under these conditions, we found 

fiP = 1.0425(25) , l/i/W = 0.50(5) , 
ft ] = 1.1075(25) , l/v^ = 0.45(5) . 

In Fig. [D] we show the overlap of the curves for the Binder cumulant obtained on 

various lattices when plotted versus (J3 - f3 c )L 1/u , with \jv fixed at 0.5 and f3 c = 1.0425. 
The quality of the overlap seems to be overall better than in the BKT scenario of Fig. [61 
but a closer inspection shows that it is worse in the region near the critical point. 

Second order transition at (3c , however, should be excluded by the numerical evidence 
that the intermediate phase is massless. In particular, two-point correlators tend to vanish 
at large distances for large volumes, whereas one should expect spontaneous symmetry 
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Figure 11: (Color online) Reduced 4-th order Binder cumulant versus {(3 — P^L 1 ^ 
for (3 C = 1.0425 and v = 2 on lattices with L ranging from 128 to 1024. 



breaking and non- vanishing values of long distance correlations in the case of second order 

(2) 

phase transition. The exclusion of the second order transition at /3c seems to be a more 
subtle problem, at least on the numerical side. One should probably simulate the system 
on larger lattices to reliably distinguish the BKT scenario from the second order one, if 
we deal with the quantities studied so far. However, one could consider more traditional 
observables to determine the order of the phase transition. The BKT phase transition is 
of infinite order. In particular, it is expected that the singular part of the free energy of 
the XY model behaves like F s ~ £~ 2 . Hence, all derivatives of the free energy are analytic 
functions of the temperature. In turn, the second derivative of the free energy shows a 
finite jump if the system undergoes a second order phase transition. We have decided, 
therefore, to compute the specific heat Cy of the Z(5) model. Figure [T2l shows the result 
of simulations for various lattice sizes. As suggested by this plot, the specific heat shows 
neither a divergence nor a finite jump. We interpret this behavior as further evidence 
in favor of the infinite order phase transition. On the theoretical side, the second order 
transition seems incompatible with the analysis of the dual transformations [151 ITE] . 

As discussed in the Introduction, in a recent work [19] it has been claimed that the 
phase transition at /3<T is not a standard BKT phase transition. The main tool in the 
analysis of Ref. [19] was the helicity modulus T, which was not considered in this work. 
The key observation in Ref. [19] was that the helicity modulus does not jump to zero 
across the phase transition. This property prompted the authors of Ref. [19J to conclude 
that the phase transition at /3c is a weaker cousin of the standard BKT transition. Our 
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Figure 12: (Color online) Specific heat Cy measured on lattices with L ranging from 16 



data indicate that there seem to be no influence of such behavior of the helicity modulus 
on other characteristic features of the phase transition. Most important on our opinion is 
the fact that both standard and Villain formulation are still in the same universality class 
since they show equal critical indices. Moreover, one could consider the critical index 
which governs the behavior of the helicity modulus [29] 




Our data for the indices (3 and rf 1 ' are compatible with a vanishing value of v. Thus, 
T = const right at the critical point. Obviously, v = also for the Villain model. The 
non-vanishing value of T at (3 < /3c seems to characterize rather the high-temperature 
phase than the massless BKT phase. Indeed, a physical interpretation given in Ref. [T9"] 
refers to lack of free vortices in the high-temperature phase of the Z(5) models. It would 
then be interesting and important to study the dynamics of the vortex-anti- vortex pairs in 
both formulations. This dynamics can indeed be different. This task, however is beyond 
the scope of the present paper. 



to 1024. 




T 




v = 2/3 — nv . 
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Appendix 



Consider general Z{5) vector Potts model with Boltzmann weight given by 

2tt Att 
Q(s) = l + 2x(p) cos—s + 2y(f3) cos— s. 

5 5 



Duality transformations read 

x d {p) = (1 + 2x{0) cos(27r/5) + 2y{fi) cos(47r/5))/(l + 2x(f3) + 2y(/3)) , 
y d (/3) = (1 + 2x(f3) cos(47r/5) + 2y(/3) cos(2tt/5))/(1 + 2x(f3) + 2y{(5)) . 

The initial couplings x(/3) and can be calculated from the Fourier transform of the 
original Boltzmann weight and in our case they are given by 

= (_i _ ^ + (_i + ^) e ^ + 2e2( 5+ ^)/(4 + 4e^ + 2e^ 5+VE ^) , (17) 
y (P) = (_i + ^5 - (1 + v^e^ + 2e^ (5+ ^ )/3 )/(4 + 4e^ + 2ei (5+ ^ )/3 ) . (18) 

It follows that 

Xd (f3) = e J(-5+V5)/3 } (19) 
VM = e~-^ + ^ . (20) 
Now, consider original and dual partition functions, 

Z(x(/3),y(P)) = C{p) Z(xM,y d (P)) ■ 

They have the same form and differ only by a smooth function C(/3). Suppose the original 
partition function is critical at /3c and (3c 2 \ These values correspond to (x^\y^) and 
(x^ 2 \y^). Let (% d \y d ) and {x d 2 \y d ) be the values of dual couplings at critical points. 
The interaction in original and dual partition functions is the same. Therefore, numerical 
values of the critical points in terms of (x, y) and (x d , yd) should be the same. However, 
we know from the solution of the self-dual equation (see [H]) that the self-dual point is 
not a critical point for the vector Potts model. This leaves the only one possibility, that 
at the critical points one must have 



and 
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Results for critical points reported in the text are /3c = 1.051 and (3^ = 1.1048. So, we 
easily find from (n"7j)-(l20|) 

x(1.051) = 0.466532 , x d (1.1048) = 0.46608 
1/(1.051) = 0.136626 , y d (1.1048) = 0.135523 



and 



z(1.1048) = 0.485097 , x d (1.051) = 0.483733 
1/(1.1048) = 0.149612 , y d (1.051) = 0.149378 . 
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